---
title: "2 - Testing H1 and H2"
author: "Arturo Bertero"
date: "2023-09-13"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
#clean
rm(list = ls())

# Libraries
library(pacman)
p_load(tidyverse, here, EGAnet, psychTools, conflicted, haven, lavaan, qgraph,
       stargazer, mgm, devtools, semPlot, psychonetrics) 


#notation
options(scipen=999)

#Conflicts
conflicts_prefer(dplyr::select)

```

# Input

```{r}
#Load data
IPBS = readRDS((here("Input", "IPBS_imputed.rds")))

#Filter smaller dataset
att = IPBS %>% 
  select(c(abort:immig)) 

```

# Processing

```{r}
# Define shortnames, longnames, and groupings 
shortnames <- names(att)
longnames <- c("Abortion","Euthanasia","Same-sex marriage",
               "Redistribution","Globalization","Immigration")

#remove labels
# att = sapply(att, haven::zap_labels)
```


## H1: Dimensionality

### EGA

```{r}
## EGA

# running EGA
att_EGA = EGA(att)
summary(att_EGA)

# Obj for dimensions
Totalgroup_comm <- list(
  "Cultural" = c(1:3),
  "Economic" = c(4:6)
)

Totalgroup_cols <- c("#966FD6B3","#008B8BB3")  

# Plot in the standard format
qgraph(att_EGA$network, labels = shortnames, nodeNames = longnames,  
         groups = Totalgroup_comm, color = Totalgroup_cols, 
         legend = FALSE, legend.cex = 0.33, vTrans = 255,
         borders = TRUE, vsize = 7.0, esize = 15, GLratio = 2,
         filetype = "jpg", layout = "spring", theme = "Borkulo",  
         filename = here("Output", "Article", "Figure_2"))

# Plot saved as tiff
qgraph(att_EGA$network, labels = shortnames, nodeNames = longnames,  
         groups = Totalgroup_comm, color = Totalgroup_cols, 
         legend = FALSE, legend.cex = 0.33, vTrans = 255,
         borders = TRUE, vsize = 7.0, esize = 15, GLratio = 2,
         filetype = "tiff", layout = "spring", theme = "Borkulo",  
         filename = here("Output", "Article", "Figure_2_tiff"))

# Boot_EGA
boot_EGA = bootEGA(att, seed = 123, iter = 1000)

#EGA compare
compare_EGA = compare.EGA.plots(att_EGA, boot_EGA, 
                                labels = "Empirical", "Bootstrap",
                                node.color = c("#966FD6B3","#008B8BB3"),
                                edge.color = c("blue", # positive edges first 
                                               "red")) 

# Stability plot
stab_EGA = dimensionStability(boot_EGA)
stab_EGA = plot(itemStability(boot_EGA))

#save
 ggsave(
      here("Output", "Supplement", "Fig_1.jpg"),
      stab_EGA,
      width = 6, height = 6, dpi = 600
    )
```

### CFA

```{r}
# Define CFA models
cfa_1_model <- 'F1 =~ abort + eutha + marria + redis + globa + immig'

cfa_2_model <- 'F1 =~ abort + eutha + marria
                F2 =~ redis + globa + immig' 

# Fit
cfa_1_fit <- cfa(cfa_1_model, data = att, std.lv = TRUE)
cfa_2_fit <- cfa(cfa_2_model, data = att, std.lv = TRUE)

# Stats
  fit_indices <- tibble(
    Model = c("One-Factor", "Two-Factor"),
    CFI = c(fitMeasures(cfa_1_fit, "cfi"), fitMeasures(cfa_2_fit, "cfi")),
    TLI = c(fitMeasures(cfa_1_fit, "tli"), fitMeasures(cfa_2_fit, "tli")),
    RMSEA = c(fitMeasures(cfa_1_fit, "rmsea"), fitMeasures(cfa_2_fit, "rmsea")),
    SRMR = c(fitMeasures(cfa_1_fit, "srmr"), fitMeasures(cfa_2_fit, "srmr"))
  )

# Test
chi2_test <- lavTestLRT(cfa_1_fit, cfa_2_fit)
chi2_test = as.data.frame(chi2_test)

# Save
stargazer(fit_indices, type = "html", summary = FALSE, rownames = FALSE, 
         out = here("Output", "Supplement", "Tab_2.doc"))

stargazer(chi2_test, type = "html", summary = FALSE, rownames = FALSE, 
         out = here("Output", "Supplement", "tab_3.doc"))
```

## H2: Network vs latent

```{r}
# CFA 2 factors

# Extract the covariance matrix
cov_att <- cov(att, use = "pairwise.complete.obs")


# Fit the CFA 2-factor model in Psychonetrics
lambda_cfa <- matrix(0, 6, 2)
lambda_cfa[1:3, 1] <- 1  # Cultural
lambda_cfa[4:6, 2] <- 1  # Economic

colnames(lambda_cfa) <- c("Cultural", "Economic")
rownames(lambda_cfa) <- colnames(att)

cfa_model <- lvm(covs = cov_att, 
                 lambda = lambda_cfa, 
                 nobs = nrow(att),
                 identification = "variance") %>%
             runmodel()

#psychonetrics::fit(cfa_model)

```

```{r}
# GGM 

# Binary adjacency matrix from EGA structure 
omega_mat <- 1 * (att_EGA$network != 0)

ggm_model <- ggm(covs = cov_att, 
                 omega = omega_mat,
                 nobs = nrow(att)) %>%
             runmodel()

```

```{r}
# Saturated model

saturated_model <- ggm(
  covs = cov(att, use = "pairwise.complete.obs"),
  omega = "full",
  nobs = nrow(att)
) %>%
  runmodel()

```

```{r}
# Compare fit (net vs latent)
compare_models <- psychonetrics::compare(cfa = cfa_model, network = ggm_model)

# Compare all three models
compare_models3 = psychonetrics::compare(
  saturated = saturated_model,
  network = ggm_model,
  cfa = cfa_model
)

# Save
stargazer(compare_models3, type = "html", summary = FALSE, rownames = FALSE, 
         out = here("Output", "Article", "Table_2.doc"))
```

## Robustness

```{r}
#Data
robustness <- att %>% 
  select(-globa)  # exclude globalization item

#EGA again

# EGA
ega_rob <- EGA(robustness)

#CFA with one dimension
lambda_cfa_rob <- matrix(1, 5, 1)  # all 5 items load on one factor
colnames(lambda_cfa_rob) <- "General"
rownames(lambda_cfa_rob) <- colnames(robustness)

cfa_rob_model <- lvm(
  covs = cov(robustness, use = "pairwise.complete.obs"),
  lambda = lambda_cfa_rob,
  nobs = nrow(robustness),
  identification = "variance"
) %>% runmodel()

#GGM
omega_rob <- 1 * (ega_rob$network != 0)

ggm_rob_model <- ggm(
  covs = cov(robustness, use = "pairwise.complete.obs"),
  omega = omega_rob,
  nobs = nrow(robustness)
) %>% runmodel()

# Compare
comp_robustness = psychonetrics::compare(cfa = cfa_rob_model, network = ggm_rob_model)

# Save
stargazer(comp_robustness, type = "html", summary = FALSE, rownames = FALSE, 
         out = here("Output", "Supplement", "Tab_4.doc"))
```

